import yt

import numpy as np

L = [-1, -1, -1]
R = [1, 1, 1]
oct_mask = np.zeros(33) # 5 refined values gives 7 * 4 + 5 octs to mask
oct_mask[[0,  5,  7, 16]] = 8
octree_mask = np.array(oct_mask, dtype=np.uint8)
quantities = {}
rng = np.random.default_rng(12345)
quantities["gas", "density"] = lambda: rng.random((29, 1))
quantities["gas", "dinos"] = (lambda: rng.random((29, 1)), "Msun")

quantities["io", "particle_position_x"] = (rng.uniform(L[0], R[0], size=100), "code_length")
quantities["io", "particle_position_y"] = (rng.uniform(L[1], R[1], size=100), "code_length")
quantities["io", "particle_position_z"] = (rng.uniform(L[2], R[2], size=100), "code_length")
quantities["io", "particle_mass"] = (rng.random(100), "Msun")
bbox = np.array([L, R]).T

ds = yt.load_octree(
    octree_mask=octree_mask,
    data=quantities,
    bbox=bbox,
    num_zones=1,
    partial_coverage=0,
)

p = yt.ProjectionPlot(ds, "x", ("gas", "density"))
p.annotate_particles(100, ptype="io")
p.save()